#include <stdio.h>


#define INDEX(x, y, z)  ((x) + (long long) ldx * ((y) + (long long) ldy * (z)))

void stencil_order_4_kern_athread(double* dst, const double* src, const double* coeffc2_dist,
                                  int start_x, int start_y, int start_z,
                                  int end_x, int end_y, int end_z, int ldx, int ldy) {
    int i, j, k;


    for(k = start_z; k < end_z; k++) {
        for(j = start_y; j < end_y; j++) {

            for(i = start_x; i < end_x; i++) {
                double coeffc2 = coeffc2_dist[INDEX(i, j, k)];

                dst[INDEX(i, j, k)] = -dst[INDEX(i, j, k)]
                                      + (2 - 7.5 * coeffc2) * src[INDEX(i, j, k)]
                                      + coeffc2 * (src[INDEX(i - 1, j, k)] * 4 / 3
                                                   + src[INDEX(i + 1, j, k)] * 4 / 3
                                                   - src[INDEX(i - 2, j, k)] / 12
                                                   - src[INDEX(i + 2, j, k)] / 12

                                                   + src[INDEX(i, j - 1, k)] * 4 / 3
                                                   + src[INDEX(i, j + 1, k)] * 4 / 3
                                                   - src[INDEX(i, j - 2, k)] / 12
                                                   - src[INDEX(i, j + 2, k)] / 12

                                                   + src[INDEX(i, j, k - 1)] * 4 / 3
                                                   + src[INDEX(i, j, k + 1)] * 4 / 3
                                                   - src[INDEX(i, j, k - 2)] / 12
                                                   - src[INDEX(i, j, k + 2)] / 12
                                                  );
            }
        }
    }
}


void stencil_order_4_pml_kern_athread(double* dst, const double* src, const double* coeffc2_dist, const double* c2dtQ_dist,
                                      int start_x, int start_y, int start_z,
                                      int end_x, int end_y, int end_z, int ldx, int ldy) {
    int i, j, k;
#ifdef __DEBUG__
    printf("entered slv func\n");
#endif

    for(k = start_z; k < end_z; k++) {
#ifdef __DEBUG__
        printf("entered slv loop k\n");
#endif
        for(j = start_y; j < end_y; j++) {

            for(i = start_x; i < end_x; i++) {
                double Qc2dt = c2dtQ_dist[INDEX(i, j, k)];
                double coeffc2 = coeffc2_dist[INDEX(i, j, k)];

                dst[INDEX(i, j, k)] = ((Qc2dt - 1) * dst[INDEX(i, j, k)]
                                       + ((2 - Qc2dt * Qc2dt) - 7.5 * coeffc2) * src[INDEX(i, j, k)]
                                       + coeffc2 * (src[INDEX(i - 1, j, k)] * 4 / 3         // line 1
                                                    + src[INDEX(i + 1, j, k)] * 4 / 3       // line 2
                                                    - src[INDEX(i - 2, j, k)] / 12          // line 3
                                                    - src[INDEX(i + 2, j, k)] / 12          // line 4

                                                    + src[INDEX(i, j - 1, k)] * 4 / 3       // line 5
                                                    + src[INDEX(i, j + 1, k)] * 4 / 3       // line 6
                                                    - src[INDEX(i, j - 2, k)] / 12          // line 7
                                                    - src[INDEX(i, j + 2, k)] / 12          // line 8

                                                    + src[INDEX(i, j, k - 1)] * 4 / 3       // line 9
                                                    + src[INDEX(i, j, k + 1)] * 4 / 3       // line 10
                                                    - src[INDEX(i, j, k - 2)] / 12          // line 11
                                                    - src[INDEX(i, j, k + 2)] / 12          // line 12
                                                   )
                                      ) / (1 + Qc2dt);
            }
        }
    }
}
